﻿using System;
using System.Drawing;
using System.Drawing.Imaging;
using MathNet.Numerics;
using MentalAlchemy.Atomics;
using MathNet.Numerics.Distributions;
using System.Collections.Generic;

namespace MentalAlchemy.Molecules.Image
{
	/// <summary>
	/// Class for implementation of image filtering operations.
	/// </summary>
	public class ImageFilter
	{
		#region - Gaussian blur. -
		public static Bitmap GaussianBlur(Bitmap image, float sigma)
		{
			var res = new Bitmap(image);

			//if (radius < Const.RADIUS_MIN) { radius = Const.RADIUS_MIN; }

			// 1. calculate gaussian means.
			float[,] table0 = ImageConversion.ConvertTo2DTableFloat(image, 0),
						table1 = ImageConversion.ConvertTo2DTableFloat(image, 1),
						table2 = ImageConversion.ConvertTo2DTableFloat(image, 2);
			float[,] means0 = CalculateGaussianMeans(table0, sigma),
						means1 = CalculateGaussianMeans(table1, sigma),
						means2 = CalculateGaussianMeans(table2, sigma);

			var rect = new Rectangle(0, 0, res.Size.Width, res.Size.Height);
			BitmapData bmpData = res.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

			int bytesSize = res.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			// 2. apply calculated values.
			int width = rect.Width,
				height = rect.Height;
			int alignByte = bmpData.Stride - width * 3;
			int count = 0;
			for (int i = 0; i < height; ++i)
			{
				for (int j = 0; j < width; ++j, count += 3)
				{
					byte temp0 = Generic.ValidateIntensityValue((int)means0[i, j]),
						temp1 = Generic.ValidateIntensityValue((int)means1[i, j]),
						temp2 = Generic.ValidateIntensityValue((int)means2[i, j]);

					bytes[count] = temp0;
					bytes[count + 1] = temp1;
					bytes[count + 2] = temp2;
				}
				count += alignByte;
			}

			System.Runtime.InteropServices.Marshal.Copy(bytes, 0, bmpData.Scan0, bytesSize);
			res.UnlockBits(bmpData);

			return res;
		}

		#region - Gaussian means. -
		/// <summary>
		/// 
		/// </summary>
		/// <param name="image">Grayscale image.</param>
		/// <param name="radius"></param>
		/// <returns>2D map of means.</returns>
		public static float[,] CalculateGaussianMeans(Bitmap image, int radius)
		{
			var table = ImageConversion.ConvertTo2DTable(image, 0);
			float[,] res = CalculateGaussianMeans(table, radius);
			return res;
		}

		public static float[,] CalculateGaussianMeans(byte[,] table, int radius)
		{
			int height = table.GetLength(0), width = table.GetLength(1);
			var tempTable = new float[height, width];
			for (int i = 0; i < height; ++i)
			{
				for (int j = 0; j < width; ++j)
				{
					tempTable[i, j] = table[i, j];
				}
			}

			var res = CalculateGaussianMeans(tempTable, radius);
			return res;
		}

		public static float[,] CalculateGaussianMeans(float[,] table, int radius)
		{
			if (table.Length <= 0) { throw new Exception("[CalculateGaussianMeans] error: Empty input 2d array."); }

			int height = table.GetLength(0), width = table.GetLength(1);
			var res = new float[height, width];

			// calculate LUT.
			var rand = new NormalDistribution(0, radius);
			var lut = new float[radius + 1];
			float lutSum = 0.0f;
			for (int i = 0; i < lut.Length; ++i)
			{
				lut[i] = (float)rand.ProbabilityDensity(i);
				if (i != 0)
				{
					lutSum += lut[i] * 8 * i;
				}
				else
				{
					lutSum += lut[i];
				}
			}
			if (lutSum > 1)
			{
				lutSum = 1.0f / lutSum;
			}

			// fill in the resulting table.
			// 1. walk through the inner part of the source table.
			//float area_1 = 1.0 / (2 * radius + 1);
			//area_1 *= area_1;
			int height_radius = height - radius,
				width_radius = width - radius;
			for (int i = radius; i < height_radius; ++i)
			{
				for (int j = radius; j < width_radius; ++j)
				{
					float sum = table[i, j] * lut[0];	// center of the neighborhood.

					// cycle through all possible radiuses starting from 1.
					for (int r = 1; r <= radius; ++r)
					{
						int up = i - r,
							down = i + r,
							left = j - r,
							right = j + r;

						float tempSum = 0.0f;
						for (int k = left; k <= right; ++k) { tempSum += table[up, k]; }
						for (int k = up + 1; k <= down; ++k) { tempSum += table[k, right]; }
						for (int k = left; k < right; ++k) { tempSum += table[down, k]; }
						for (int k = up + 1; k < down; ++k) { tempSum += table[k, left]; }

						sum += tempSum * lut[r];
					}

					res[i, j] = sum * lutSum;
				}
			}
			// 2. Walk through the borders of the source table.
			// 2.1. Upper horizontal border.
			for (int i = 0; i < radius; ++i)
			{
				for (int j = 0; j < width; ++j)
				{
					int up = i - radius,
						down = i + radius,
						left = j - radius,
						right = j + radius;

					if (up < 0) up = 0;
					if (down >= height) down = height - 1;
					if (left < 0) left = 0;
					if (right >= width) right = width - 1;

					float sum = 0.0f;
					for (int i1 = up; i1 <= down; ++i1)
					{
						for (int j1 = left; j1 <= right; ++j1)
						{
							int r1 = Math.Abs(i1 - i), r2 = Math.Abs(j1 - j);
							int dist = r1 > r2 ? r1 : r2;
							sum += table[i1, j1] * lut[dist];
						}
					}
					res[i, j] = sum * lutSum;
				}
			}
			// 2.2. Rightmost vertical border.
			for (int i = radius; i < height; ++i)
			{
				for (int j = width_radius; j < width; ++j)
				{
					int up = i - radius,
						down = i + radius,
						left = j - radius,
						right = j + radius;

					if (up < 0) up = 0;
					if (down >= height) down = height - 1;
					if (left < 0) left = 0;
					if (right >= width) right = width - 1;

					float sum = 0.0f;
					for (int i1 = up; i1 <= down; ++i1)
					{
						for (int j1 = left; j1 <= right; ++j1)
						{
							int r1 = Math.Abs(i1 - i), r2 = Math.Abs(j1 - j);
							int dist = r1 > r2 ? r1 : r2;
							sum += table[i1, j1] * lut[dist];
						}
					}
					res[i, j] = sum * lutSum;
				}
			}
			// 2.3. Bottom horizontal border.
			for (int i = height_radius; i < height; ++i)
			{
				for (int j = 0; j < width_radius; ++j)
				{
					int up = i - radius,
						down = i + radius,
						left = j - radius,
						right = j + radius;

					if (up < 0) up = 0;
					if (down >= height) down = height - 1;
					if (left < 0) left = 0;
					if (right >= width) right = width - 1;

					float sum = 0.0f;
					for (int i1 = up; i1 <= down; ++i1)
					{
						for (int j1 = left; j1 <= right; ++j1)
						{
							int r1 = Math.Abs(i1 - i), r2 = Math.Abs(j1 - j);
							int dist = r1 > r2 ? r1 : r2;
							sum += table[i1, j1] * lut[dist];
						}
					}
					res[i, j] = sum * lutSum;
				}
			}
			// 2.4. Leftmost vertical border.
			for (int i = radius; i < height_radius; ++i)
			{
				for (int j = 0; j < radius; ++j)
				{
					int up = i - radius,
						down = i + radius,
						left = j - radius,
						right = j + radius;

					if (up < 0) up = 0;
					if (down >= height) down = height - 1;
					if (left < 0) left = 0;
					if (right >= width) right = width - 1;

					float sum = 0.0f;
					for (int i1 = up; i1 <= down; ++i1)
					{
						for (int j1 = left; j1 <= right; ++j1)
						{
							int r1 = Math.Abs(i1 - i), r2 = Math.Abs(j1 - j);
							int dist = r1 > r2 ? r1 : r2;
							sum += table[i1, j1] * lut[dist];
						}
					}
					res[i, j] = sum * lutSum;
				}
			}

			return res;
		}

		/// <summary>
		/// [atomic]
		/// 
		/// Non OpenMP version.
		/// </summary>
		/// <param name="table"></param>
		/// <param name="sigma"></param>
		/// <returns></returns>
		public static float[,] CalculateGaussianMeansEx(float[,] table, float sigma)
		{
			if (table.Length <= 0) { throw new Exception("[CalculateGaussianMeans] error: Empty input 2d array."); }

			var kernel = CalculateGaussianKernelFloat(sigma);
			var sum = MatrixMath.Sum(kernel);
			//int radius = kernel.GetLength(0)/2;
			kernel = MatrixMath.Mul(kernel, 1.0f / sum);
			var res = Generic.Convolve(table, kernel);
			return res;
		}

		/// <summary>
		/// [atomic]
		/// 
		/// Calculates Gaussian means for the given 2d table and sigma.
		/// Uses OpenMP.
		/// </summary>
		/// <param name="table"></param>
		/// <param name="sigma"></param>
		/// <returns></returns>
		public static float[,] CalculateGaussianMeans(float[,] table, float sigma)
		{
			if (table.Length <= 0) { throw new Exception("[CalculateGaussianMeans] error: Empty input 2d array."); }

			var kernel = CalculateGaussianKernelFloat(sigma);
			var sum = MatrixMath.Sum(kernel);
			kernel = MatrixMath.Mul(kernel, 1.0f / sum);
			var res = Generic.Convolve(table, kernel);
			return res;
		}

		/// <summary>
		/// [atomic]
		/// 
		/// Calculate Gaussian kernel for the given sigma value.
		/// </summary>
		/// <param name="sigma">Gaussian distribution sigma parameter.</param>
		/// <returns>Gaussian kernel.</returns>
		public static float[,] CalculateGaussianKernelFloat(float sigma)
		{
			// 1. define kernel size.
			var tester = new List<float>();
			const float eps = 0.001f;
			float r = 1.0f;
			var rand = new NormalDistribution(0, sigma);
			do
			{
				tester.Add((float)rand.ProbabilityDensity(r));
				++r;
			} while (Math.Abs(tester[tester.Count - 1]) > eps);

			int radius = tester.Count; int diameter = radius * 2 + 1;
			var kernel = new float[diameter, diameter];
			kernel[radius, radius] = (float)rand.ProbabilityDensity(0);
			for (int i = 0; i < radius; ++i)
			{
				for (int j = i; j < radius; ++j)
				{
					float x = radius - i, y = radius - j;
					var dist = (float)Math.Sqrt(x * x + y * y);
					var val = (float)rand.ProbabilityDensity(dist);

					kernel[i, j] = val;
					kernel[j, i] = val;

					kernel[i, diameter - j - 1] = val;
					kernel[diameter - j - 1, i] = val;

					kernel[j, diameter - i - 1] = val;
					kernel[diameter - i - 1, j] = val;

					kernel[diameter - i - 1, diameter - j - 1] = val;
					kernel[diameter - j - 1, diameter - i - 1] = val;
				}

				var val2 = (float)rand.ProbabilityDensity(radius - i);
				kernel[radius, i] = val2;
				kernel[i, radius] = val2;
				kernel[diameter - i - 1, radius] = val2;
				kernel[radius, diameter - i - 1] = val2;
			}

			return kernel;
		}
		#endregion
		#endregion

		#region - Difference Of Gaussians. -
		public class DifferenceOfGaussians
		{
			/// <summary>
			/// [molecule]
			/// 
			/// Calculates DoG using s1 / s2 ratio equal to 1.6 (by Marr and Hildrith).
			/// </summary>
			/// <param name="img"></param>
			/// <param name="radius"></param>
			/// <returns></returns>
			public static byte[,] Calculate(Bitmap img, int radius)
			{
				var temp = new Bitmap(img);
				byte[,] grey = ImageConversion.ConvertToGrayscaleTable(temp);

				var gauss1 = CalculateGaussianMeans(grey, radius);
				var gauss2 = CalculateGaussianMeans(grey, (int)(radius * 1.6 + 0.5));

				int width = img.Width, height = img.Height;
				var res = new byte[height, width];
				for (int i = 0; i < height; ++i)
				{
					for (int j = 0; j < width; ++j)
					{
						var value = (int)(gauss1[i, j] - gauss2[i, j]);
						res[i, j] = Generic.ValidateIntensityValue(value);
					}
				}
				return res;
			}
		} 
		#endregion

		#region - Laplacian. -
		// todo: write correct laplacian via calculation of laplacian kernel and convolution with this kernel.
		public static Bitmap Laplacian(Bitmap image)
		{
			byte[,] table = ImageConversion.ConvertToGrayscaleTable(image);

			// todo: calculate kernel.
			int[,] kernel = {{-1, -1, -1, -1, -1},
							{-1, -1, 1, -1, -1},
							{-1, 1, 16, 1, -1},
							{-1, -1, 1, -1, -1},
							{-1, -1, -1, -1, -1}};

			byte[,] bytes = Generic.Convolve(table, kernel);
			return ImageConversion.CreateBmpFromGrayscaleTable(bytes);
		} 
		#endregion

		#region - Laplacian Of Gaussian. -
		public class LaplacianOfGaussian
		{
			/// <summary>
			/// Calculates LoG using sigma = 1.4.
			/// </summary>
			/// <param name="img"></param>
			/// <returns></returns>
			public static byte[,] Calculate(Bitmap img)
			{
				return Calculate(img, 1.4f);
			}

			/// <summary>
			/// Calculates LoG using arbitrary sigma.
			/// </summary>
			/// <param name="img"></param>
			/// <param name="sigma"></param>
			/// <returns></returns>
			public static byte[,] Calculate(Bitmap img, float sigma)
			{
				if (sigma <= 0.0) { throw new Exception("[LaplacianOfGaussian :: Calculate] exception: Wrong sigma value."); }

				// calculate convolution kernel.
				var kernel = CalculateKernelFloat(sigma);

				var gray = ImageConversion.ConvertToGrayscaleTable(img);
				var log = Generic.Convolve(gray, kernel);
				return log;
			}

			/// <summary>
			/// Calculates LoG using arbitrary sigma.
			/// The resulting bytes array has brightness offset equal to 128 for display purposes.
			/// </summary>
			/// <param name="img"></param>
			/// <param name="sigma"></param>
			/// <returns></returns>
			public static byte[,] CalculateOffset(Bitmap img, float sigma)
			{
				if (sigma <= 0.0) { throw new Exception("[LaplacianOfGaussian :: Calculate] exception: Wrong sigma value."); }

				// calculate convolution kernel.
				var kernel = CalculateKernelFloat(sigma);

				var gray = ImageConversion.ConvertToGrayscaleTable(img);
				var log = Generic.ConvolveOffset(gray, kernel, 128);
				return log;
			}

			/// <summary>
			/// [molecule]
			/// 
			/// Calculates LoG using arbitrary sigma and operating with floating point values (instead of bytes).
			/// The resulting array has no offset since floats allow negative values.
			/// </summary>
			/// <param name="img">Input image.</param>
			/// <param name="sigma">Sigma for LoG.</param>
			/// <param name="offset">Offset, which is added to each resulting value.</param>
			/// <returns>Resulting LoG 2d map.</returns>
			public static float[,] CalculateOffsetFloat(Bitmap img, float sigma, byte offset)
			{
				if (sigma <= 0.0) { throw new Exception("[LaplacianOfGaussian :: Calculate] exception: Wrong sigma value."); }

				// calculate convolution kernel.
				var kernel = CalculateKernelFloat(sigma);

				var gray = ImageConversion.ConvertToGrayscaleTableFloat(img);
				var log = Generic.ConvolveOffset(gray, kernel, offset);
				return log;
			}

			public static float[,] CalculateKernelFloat(float sigma)
			{
				// 1. define kernel size.
				var tester = new List<float>();
				//float normC = 1.0f,
				const float eps = 0.001f;
				float x = 1.0f;
				do
				{
					tester.Add(LoGResponce(x, 0, sigma));
					++x;
				} while (Math.Abs(tester[tester.Count - 1]) > eps);

				int radius = tester.Count; int diameter = radius * 2 + 1;
				var kernel = new float[diameter, diameter];
				kernel[radius, radius] = LoGResponce(0, 0, sigma);
				for (int i = 0; i < radius; ++i)
				{
					for (int j = i; j < radius; ++j)
					{
						float val = LoGResponce(radius - i, radius - j, sigma);

						kernel[i, j] = val;
						kernel[j, i] = val;

						kernel[i, diameter - j - 1] = val;
						kernel[diameter - j - 1, i] = val;

						kernel[j, diameter - i - 1] = val;
						kernel[diameter - i - 1, j] = val;

						kernel[diameter - i - 1, diameter - j - 1] = val;
						kernel[diameter - j - 1, diameter - i - 1] = val;
					}

					float val2 = LoGResponce(radius - i, 0, sigma);
					kernel[radius, i] = val2;
					kernel[i, radius] = val2;
					kernel[diameter - i - 1, radius] = val2;
					kernel[radius, diameter - i - 1] = val2;
				}

				float sum = MatrixMath.Sum(kernel);
				kernel[radius, radius] -= sum;
				//sum = EvoMath.MatrixMath.Sum(kernel);
				//kernel = EvoMath.Mul(kernel, normC);

				return kernel;
			}

			public static float LoGResponce(float x, float y, float sigma)
			{
				float sigma2 = sigma * sigma, sigma2_1 = 1 / sigma2;
				float r2 = x * x + y * y;
				float pow = 0.5f * r2 * sigma2_1;
				float res = -sigma2_1 * sigma2_1 / (float)Math.PI;
				res *= (1 - pow);
				res *= (float)Math.Exp(-pow);
				return res;
			}
		} 
		#endregion

		#region - Edge Detector class. -
		public class EdgeDetector
		{
			public static int Sobel(Bitmap image, float argT)
			{
				switch (image.PixelFormat)
				{
					case PixelFormat.Format24bppRgb:
						Sobel24bit(image, argT);
						break;
					case PixelFormat.Format32bppArgb:
						Sobel24bit(image, argT);
						break;
					case PixelFormat.Format8bppIndexed:
						Sobel8bit(image, argT);
						break;
					default:
						throw new Exception("EdgeDetector::Sobel: Unsupported pixel format (" + image.PixelFormat + ").");
				}
				return 0;
			}

			public static int Sobel24bit(Bitmap image, float argT)
			{
				if (image != null)
				{
					int width = image.Width, height = image.Height;
					var d = new float[2];

					var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
					BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
					//int bytesSize = image.Height * image.Width * 3;
					int bytesSize = image.Height * bmpData.Stride;
					var bytes = new byte[bytesSize];
					System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

					var newBytes = new byte[bytesSize];

					// Apply Sobel's transformation (see K.Lindley for details)
					// a b c
					// d e f
					// g h i
					int alignByte = bmpData.Stride - width * 3;
					int stride = bmpData.Stride;
					int count = bmpData.Stride;

					for (int i = 1; i < (height - 1); ++i)
					{
						newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						count += 3;	// skip the first pixel
						//int curRow = i * width * 3 + alignByte,
						//    prevRow = (i - 1) * width * 3 + alignByte,
						//    nextRow = (i + 1) * width * 3 + alignByte;
						int curRow = i * stride,
							prevRow = curRow - stride,
							nextRow = curRow + stride;
						for (int j = 1; j < (width - 1); ++j, count += 3)
						{
							int curCol = j * 3,
								prevCol = curCol - 3,
								nextCol = curCol + 3;

							var av1 = bytes[prevRow + curCol] + bytes[prevRow + nextCol] + bytes[curRow + nextCol];	// sum of b-c-f ok
							var av2 = bytes[nextRow + curCol] + bytes[nextRow + prevCol] + bytes[curRow + prevCol];	// sum of d-g-h ok
							d[0] = Math.Abs(av1 - av2);	// for a-e-i

							av1 = bytes[prevRow + prevCol] + bytes[curRow + prevCol] + bytes[nextRow + prevCol];	// sum of a-d-g
							av2 = bytes[nextRow + nextCol] + bytes[curRow + nextCol] + bytes[prevRow + nextCol];	// sum of c-f-i
							d[1] = Math.Abs(av1 - av2);	// for b-e-h
							if (d[1] > d[0])
							{
								d[0] = d[1];
							}

							av1 = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[curRow + prevCol];	// sum of a-b-d ok
							av2 = bytes[nextRow + nextCol] + bytes[nextRow + curCol] + bytes[curRow + nextCol];	// sum of f-h-i ok
							d[1] = Math.Abs(av1 - av2);	// for c-e-g
							if (d[1] > d[0])
							{
								d[0] = d[1];
							}

							av1 = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[prevRow + nextCol];	// sum of a-b-c
							av2 = bytes[nextRow + nextCol] + bytes[nextRow + curCol] + bytes[nextRow + prevCol];	// sum of g-h-i
							d[1] = Math.Abs(av1 - av2);	// for d-e-f
							if (d[1] > d[0])
							{
								d[0] = d[1];
							}

							d[0] = d[0] / 3.0f + 0.5f;	// normalization and preparation for proper rounding
							byte temp;
							if (d[0] < argT)
							{
								temp = 0;
							}
							else
							{
								temp = (byte)d[0];
							}
							newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = temp;
							//image.SetPixel(j, i, Color.FromArgb(temp, temp, temp));
						}
						newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						count += 3 + alignByte;	// skip the last pixel

						//if (width % 2 != 0)
						//{
						//    //newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						//    ++count;	// skip one more pixel for alignment.
						//}
					}

					System.Runtime.InteropServices.Marshal.Copy(newBytes, 0, bmpData.Scan0, bytesSize);
					image.UnlockBits(bmpData);

					return 1;

				}
				//cout<<"edgeSobel:\tCan't transform empty image\n";
				return 0;
			}

			public static int Sobel8bit(Bitmap image, float argT)
			{
				if (image != null)
				{
					int width = image.Width, height = image.Height;
					var d = new float[2];

					var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
					BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
					int bytesSize = image.Height * bmpData.Stride;
					var bytes = new byte[bytesSize];
					System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

					var newBytes = new byte[bytesSize];

					// Apply Sobel's transformation (see K.Lindley for details)
					// a b c
					// d e f
					// g h i
					int alignByte = bmpData.Stride - width;
					int stride = bmpData.Stride;
					int count = width + alignByte;
					for (int i = 1; i < (height - 1); ++i)
					{
						newBytes[count] = 0;
						++count;	// skip the first pixel
						int curRow = i * stride,
							prevRow = (i - 1) * stride,
							nextRow = (i + 1) * stride;
						//int curRow = i * width,
						//    prevRow = (i - 1) * width,
						//    nextRow = (i + 1) * width;
						for (int j = 1; j < (width - 1); ++j, ++count)
						{
							int curCol = j,
								prevCol = (j - 1),
								nextCol = (j + 1);

							var av1 = bytes[prevRow + curCol] + bytes[prevRow + nextCol] + bytes[curRow + nextCol];	// sum of b-c-f ok
							var av2 = bytes[nextRow + curCol] + bytes[nextRow + prevCol] + bytes[curRow + prevCol];	// sum of d-g-h ok
							d[0] = Math.Abs(av1 - av2);	// for a-e-i

							av1 = bytes[prevRow + prevCol] + bytes[curRow + prevCol] + bytes[nextRow + prevCol];	// sum of a-d-g
							av2 = bytes[nextRow + nextCol] + bytes[curRow + nextCol] + bytes[prevRow + nextCol];	// sum of c-f-i
							d[1] = Math.Abs(av1 - av2);	// for b-e-h
							if (d[1] > d[0])
							{
								d[0] = d[1];
							}

							av1 = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[curRow + prevCol];	// sum of a-b-d ok
							av2 = bytes[nextRow + nextCol] + bytes[nextRow + curCol] + bytes[curRow + nextCol];	// sum of f-h-i ok
							d[1] = Math.Abs(av1 - av2);	// for c-e-g
							if (d[1] > d[0])
							{
								d[0] = d[1];
							}

							av1 = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[prevRow + nextCol];	// sum of a-b-c
							av2 = bytes[nextRow + nextCol] + bytes[nextRow + curCol] + bytes[nextRow + prevCol];	// sum of g-h-i
							d[1] = Math.Abs(av1 - av2);	// for d-e-f
							if (d[1] > d[0])
							{
								d[0] = d[1];
							}

							//d[0] = d[0] / 3.0 + 0.5;	// normalization and preparation for proper rounding
							byte temp;
							if (d[0] < argT)
							{
								temp = 0;
							}
							else
							{
								temp = (byte)d[0];
							}
							newBytes[count] = temp;
							//image.SetPixel(j, i, Color.FromArgb(temp, temp, temp));
						}
						newBytes[count] = 0;
						count += 1 + alignByte; // skip the last pixel
					}

					System.Runtime.InteropServices.Marshal.Copy(newBytes, 0, bmpData.Scan0, bytesSize);
					image.UnlockBits(bmpData);

					return 1;
				}

				return 0;
			}

			public static int SobelMod(Bitmap image, ref Struct.EdgeInfo edgeInfo)
			{
				if (image != null)
				{
					int width = image.Width, height = image.Height;
					var d = new float[2];

					var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
					BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
					int bytesSize = image.Height * bmpData.Stride;
					var bytes = new byte[bytesSize];
					System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

					var newBytes = new byte[bytesSize];
					var edgeFreq = new int[Const.GRAYLEVELS_COUNT];
					int weakEdgeCount = 0, mediumEdgeCount = 0, strongEdgeCount = 0;

					// reset edge info structure.
					edgeInfo.EdgeFreq = new int[Const.GRAYLEVELS_COUNT];
					edgeInfo.WeakEdgeCount = edgeInfo.MediumEdgeCount = edgeInfo.StrongEdgeCount = 0;

					// Apply Sobel's transformation (see K.Lindley for details)
					// a b c
					// d e f
					// g h i
					int alignByte = bmpData.Stride - width * 3;
					int stride = bmpData.Stride;
					int count = bmpData.Stride; // proceed to the 2-nd row.
					for (int i = 1; i < (height - 1); ++i)
					{
						newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						count += 3;	// skip the first pixel
						int curRow = i * stride,
							prevRow = curRow - stride,
							nextRow = curRow + stride;
						for (int j = 1; j < (width - 1); ++j, count += 3)
						{
							int curCol = j * 3,
								prevCol = curCol - 3,
								nextCol = curCol + 3;

							// calculate sum of neighboring pixels.
							int nSum = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[prevRow + nextCol]
										+ bytes[curRow + prevCol] + bytes[curRow + nextCol]
										+ bytes[nextRow + prevCol] + bytes[nextRow + curCol] + bytes[nextRow + nextCol];
							// calculate brightness rate.
							nSum /= 8;
							//float bRate = (bytes[curRow + curCol] > nSum)?
							//                (float)bytes[curRow + curCol] / (float)nSum
							//                :
							//                (float)nSum / (float)bytes[curRow + curCol];
							//float bRate = (float)Math.Abs(bytes[curRow + curCol] - nSum) / (float)nSum;
							float bRate = (bytes[curRow + curCol] - nSum) / (float)nSum;
							if (bRate < 0) { bRate = -bRate; }

							byte temp = 0;

							float aspect = 1.0f;
							if (nSum < 97.0f) { aspect = 97.0f / nSum; }
							else if (nSum > 159.0f) { aspect = nSum / 157.0f; }
							if (bRate < Const.BRIGHTNESS_RATE * aspect)
							{
								//edgeInfo.WeakEdgeCount++;
								weakEdgeCount++;
							}
							else
							{
								var av1 = bytes[prevRow + curCol] + bytes[prevRow + nextCol] + bytes[curRow + nextCol];	// sum of b-c-f ok
								var av2 = bytes[nextRow + curCol] + bytes[nextRow + prevCol] + bytes[curRow + prevCol];	// sum of d-g-h ok
								//d[0] = Math.Abs(av1 - av2);	// for a-e-i
								d[0] = av1 - av2;	// for a-e-i
								if (d[0] < 0) { d[0] = -d[0]; }

								av1 = bytes[prevRow + prevCol] + bytes[curRow + prevCol] + bytes[nextRow + prevCol];	// sum of a-d-g
								av2 = bytes[nextRow + nextCol] + bytes[curRow + nextCol] + bytes[prevRow + nextCol];	// sum of c-f-i
								//d[1] = Math.Abs(av1 - av2);	// for b-e-h
								d[1] = av1 - av2;	// for b-e-h
								if (d[1] < 0) { d[1] = -d[1]; }
								if (d[1] > d[0])
								{
									d[0] = d[1];
								}

								av1 = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[curRow + prevCol];	// sum of a-b-d ok
								av2 = bytes[nextRow + nextCol] + bytes[nextRow + curCol] + bytes[curRow + nextCol];	// sum of f-h-i ok
								//d[1] = Math.Abs(av1 - av2);	// for c-e-g
								d[1] = av1 - av2;	// for b-e-h
								if (d[1] < 0) { d[1] = -d[1]; }
								if (d[1] > d[0])
								{
									d[0] = d[1];
								}

								av1 = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[prevRow + nextCol];	// sum of a-b-c
								av2 = bytes[nextRow + nextCol] + bytes[nextRow + curCol] + bytes[nextRow + prevCol];	// sum of g-h-i
								//d[1] = Math.Abs(av1 - av2);	// for d-e-f
								d[1] = av1 - av2;	// for b-e-h
								if (d[1] < 0) { d[1] = -d[1]; }
								if (d[1] > d[0])
								{
									d[0] = d[1];
								}

								d[0] = d[0] / 3.0f + 0.5f;	// normalization and preparation for proper rounding

								if (bRate < Const.BRIGHTNESS_RATE_EXPMUL)
								{
									temp = (byte)d[0];
									//edgeInfo.MediumEdgeCount++;
									mediumEdgeCount++;
								}
								else
								{
									//int sqr = (int)(d[0] * d[0]);
									//temp = (sqr <= Const.GRAYLEVEL_MAX) ? (byte)sqr : (byte)Const.GRAYLEVEL_MAX;
									temp = (byte)d[0];
									//edgeInfo.StrongEdgeCount++;
									strongEdgeCount++;
								}
							}

							//edgeInfo.EdgeFreq[temp]++;
							edgeFreq[temp]++;
							newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = temp;
							//image.SetPixel(j, i, Color.FromArgb(temp, temp, temp));
						}
						newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						count += 3 + alignByte;	// skip the last pixel
					}	// for (i = 1; i < (height - 1); ++i)

					edgeInfo.EdgeFreq = edgeFreq;
					edgeInfo.WeakEdgeCount = weakEdgeCount;
					edgeInfo.MediumEdgeCount = mediumEdgeCount;
					edgeInfo.StrongEdgeCount = strongEdgeCount;

					System.Runtime.InteropServices.Marshal.Copy(newBytes, 0, bmpData.Scan0, bytesSize);
					image.UnlockBits(bmpData);

					float area_1 = 1.0f / (width * height);
					edgeInfo.WeakEdgeRate = edgeInfo.WeakEdgeCount * area_1;
					edgeInfo.MediumEdgeRate = edgeInfo.MediumEdgeCount * area_1;
					edgeInfo.StrongEdgeRate = edgeInfo.StrongEdgeCount * area_1;

					return 1;
				}

				return 0;
			}

			public static int Phaeco(Bitmap image, float argT)
			{
				if (image != null)
				{
					int width = image.Width, height = image.Height;

					var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
					BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
					//int bytesSize = image.Height * image.Width * 3;
					int bytesSize = image.Height * bmpData.Stride;
					var bytes = new byte[bytesSize];
					System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

					var newBytes = new byte[bytesSize];

					// Apply Sobel's transformation (see K.Lindley for details)
					// a b c
					// d e f
					// g h i
					int alignByte = bmpData.Stride - width * 3;
					int count = bmpData.Stride; // proceed to the 2-nd row.

					for (int i = 1; i < (height - 1); ++i)
					{
						newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						count += 3;	// skip the first pixel
						int curRow = i * bmpData.Stride,
							prevRow = curRow - bmpData.Stride,
							nextRow = curRow + bmpData.Stride;
						for (int j = 1; j < (width - 1); ++j, count += 3)
						{
							int curCol = j * 3,
								prevCol = curCol - 3,
								nextCol = curCol + 3;

							int sum = bytes[prevRow + prevCol] + bytes[prevRow + curCol] + bytes[prevRow + nextCol]
									+ bytes[curRow + prevCol] - 8 * bytes[curRow + curCol] + bytes[curRow + nextCol]
									+ bytes[nextRow + prevCol] + bytes[nextRow + curCol] + bytes[nextRow + nextCol];
							sum = -sum;
							if (sum < Const.GRAYLEVEL_MIN)
							{
								sum = Const.GRAYLEVEL_MIN;
							}
							else if (sum > Const.GRAYLEVEL_MAX)
							{
								sum = Const.GRAYLEVEL_MAX;
							}
							newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = (byte)(sum);
						}
						newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						count += 3 + alignByte;	// skip the last pixel

						//if (width % 2 != 0)
						//{
						//    //newBytes[count] = newBytes[count + 1] = newBytes[count + 2] = 0;
						//    ++count;	// skip one more pixel for alignment.
						//}
					}

					System.Runtime.InteropServices.Marshal.Copy(newBytes, 0, bmpData.Scan0, bytesSize);
					image.UnlockBits(bmpData);

					return 1;

				}

				return 0;
			}
		}
		#endregion

		#region - 2D Gabor Filter class. -
		/// <summary>
		/// [molecule]
		/// 
		/// Class to implement calculation of the Gabor filter mask of
		/// prescribed size and parameters. See 'Tutorial on Gabor Filters' by
		/// Javier R. Movellan for fomulae.
		/// </summary>
		public class GaborFilter
		{
			/// <summary>
			/// [molecule]
			/// 
			/// Creates mask of the given size for Gabor filter with the prescribed
			/// parameters.
			/// </summary>
			/// <param name="maskSize">Size of the mask.</param>
			/// <param name="hFreq">Horizontal frequency.</param>
			/// <param name="vFreq">Vertical frequency.</param>
			/// <param name="phase">Phase shift.</param>
			/// <param name="geParams">Parameters for the 2D Gaussian envelope.</param>
			/// <returns>Calculated mask array.</returns>
			public static float[,] CreateFilterMask (Size maskSize, float hFreq, float vFreq, float phase, GaussianEnvelopeParameters geParams)
			{
				if (maskSize.Height % 2 == 0 || maskSize.Width % 2 == 0)
				{
					throw new Exception("[GaborFilter.CreateFilterMask]: Invalid mask size: dimensions should be odd.");
				}

				var res = new float[maskSize.Height, maskSize.Width];
				var horRadius = maskSize.Width / 2;
				var verRadius = maskSize.Height / 2;

				for (int i = 0; i < maskSize.Height; ++i )
				{
					float y = verRadius - i;	// flip Oy
					for (int j = 0; j < maskSize.Width; j++)
					{
						float x = j - horRadius;
						var point = new PointF(x, y);

						// 1. calculate carrier.
						var carrier = CalculateCarrier(point, hFreq, vFreq, phase);

						// 2. calculate envelope.
						var env = CalculateEnvelope(point, geParams);

						res[i, j] = (float)(carrier.Real*env);
					}
				}

				return res;
			}

			/// <summary>
			/// [molecule]
			/// 
			/// Creates scale-free mask of the given size for Gabor filter with the prescribed
			/// parameters. The mask is formed via sampling from square with side = 2 and which
			/// center is located at origin. As a result masks of different sizes but with
			/// the same parameters will just represent the same mask but at different spatial
			/// resolutions.
			/// </summary>
			/// <param name="maskSize">Size of the mask.</param>
			/// <param name="hFreq">Horizontal frequency.</param>
			/// <param name="vFreq">Vertical frequency.</param>
			/// <param name="phase">Phase shift.</param>
			/// <param name="geParams">Parameters for 2D Gaussian envelope.</param>
			/// <returns>Calculated mask array.</returns>
			public static float[,] CreateFilterMaskScaleFree(Size maskSize, float hFreq, float vFreq, float phase, GaussianEnvelopeParameters geParams)
			{
				if (maskSize.Height % 2 == 0 || maskSize.Width % 2 == 0)
				{
					throw new Exception("[GaborFilter.CreateFilterMask]: Invalid mask size: dimensions should be odd.");
				}
				if (maskSize.Height <= 1 || maskSize.Width <= 1)
				{
					throw new Exception("[GaborFilter.CreateFilterMask]: Invalid mask size: The mask is too small.");
				}

				var res = new float[maskSize.Height, maskSize.Width];
				float xStep = 2.0f/(maskSize.Width - 1), xStep_1 = 1.0f / xStep;
				float yStep = 2.0f / (maskSize.Height - 1), yStep_1 = 1.0f / yStep;

				for (float y = -1; y <= 1; y += yStep)
				{
					float i = maskSize.Height - (y + 1) * yStep_1 - 1;	// flip Oy
					for (float x = -1; x <= 1 ; x += xStep)
					{
						float j = (x + 1) * xStep_1;
						var point = new PointF(x, y);

						// 1. calculate carrier.
						Complex carrier = CalculateCarrier(point, hFreq, vFreq, phase);

						// 2. calculate envelope.
						float env = CalculateEnvelope(point, geParams);

						res[(int)(i + 0.5f), (int)(j + 0.5f)] = (float)(carrier.Real * env);
					}
				}

				return res;
			}

			/// <summary>
			/// [molecule]
			/// 
			/// Calculate carrier value at the given point with specified parameters.
			/// </summary>
			/// <param name="point">Point coordinates.</param>
			/// <param name="hFreq">Horizontal frequency.</param>
			/// <param name="vFreq">Vertical frequency.</param>
			/// <param name="phase">Phase shift.</param>
			/// <returns>Carrier value.</returns>
			public static Complex CalculateCarrier (PointF point, float hFreq, float vFreq, float phase)
			{
				var res = new Complex();

				var arg = (float)(2 * Math.PI * (hFreq*point.X + vFreq*point.Y) + phase);
				res.Real = Math.Cos(arg);
				res.Imag = Math.Sin(arg);
				return res;
			}

			/// <summary>
			/// Calculate envelope value for Gabor filter mask calculation.
			/// </summary>
			/// <param name="point">Point coordinates to calculate envelope value.</param>
			/// <param name="geParams">Parameters of Gaussian envelope.</param>
			/// <returns>Resulting envelope value.</returns>
			public static float CalculateEnvelope (PointF point, GaussianEnvelopeParameters geParams)
			{
				var sinT = (float)Math.Sin(geParams.Angle);
				var cosT = (float)Math.Cos(geParams.Angle);

				float dx = point.X - geParams.Center.X,
						dy = point.Y - geParams.Center.Y;

				float xRot = dx*cosT + dy*sinT,
				      yRot = -dx*sinT + dy*cosT;

				float arg = geParams.ScaleXSqr * xRot * xRot
							+ geParams.ScaleYSqr * yRot * yRot;
				var res = (float)(geParams.Magnitude*Math.Exp(-Math.PI * (arg)));

				return res;
			}
		}
		#endregion

		#region - Max. -
		/// <summary>
		/// Leaves table value only if neighboring cells are all having lesser value.
		/// </summary>
		/// <param name="table"></param>
		/// <returns></returns>
		public static float[,] Max(float[,] table)
		{
			int height = table.GetLength(0), width = table.GetLength(1);
			var res = new float[height, width];
			for (int i=1; i<height-1; ++i)
			{
				for (int j = 1; j < width-1; ++j)
				{
					float temp = table[i, j];
					// previous row.
					if (temp <= table[i - 1, j - 1])
					{
						res[i, j] = 0.0f;
						continue;
					}
					if (temp <= table[i - 1, j])
					{
						res[i, j] = 0.0f;
						continue;
					}
					if (temp <= table[i - 1, j + 1])
					{
						res[i, j] = 0.0f;
						continue;
					}
					// current row.
					if (temp <= table[i, j - 1])
					{
						res[i, j] = 0.0f;
						continue;
					}
					if (temp <= table[i, j + 1])
					{
						res[i, j] = 0.0f;
						continue;
					}
					// next row.
					if (temp <= table[i + 1, j - 1])
					{
						res[i, j] = 0.0f;
						continue;
					}
					if (temp <= table[i + 1, j])
					{
						res[i, j] = 0.0f;
						continue;
					}
					if (temp <= table[i + 1, j + 1])
					{
						res[i, j] = 0.0f;
						continue;
					}

					res[i, j] = temp;
				}
			}

			return res;
		} 
		#endregion

		#region - Binarization. -
		public static Bitmap Binarize(Bitmap src, float argT)
		{
			var res = new Bitmap(src);

			var rect = new Rectangle(0, 0, res.Size.Width, res.Size.Height);
			BitmapData bmpData = res.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

			int bytesSize = res.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			for (int i = 0; i < bytesSize; ++i)
			{
				bytes[i] = (bytes[i] > argT) ? (byte)Const.GRAYLEVEL_MAX : (byte)Const.GRAYLEVEL_MIN;
			}

			System.Runtime.InteropServices.Marshal.Copy(bytes, 0, bmpData.Scan0, bytesSize);
			res.UnlockBits(bmpData);

			return res;
		}
		#endregion

		#region - Color filtering. -
		/// <summary>
		/// Leave only those colors which are listed in a given array.
		/// </summary>
		/// <param name="image"></param>
		/// <param name="colors"></param>
		public static void LeaveColors(Bitmap image, List<Color> colors)
		{
			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
			int bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			for (int i = 0; i < bytesSize; i += 3)
			{
				if (!colors.Contains(Color.FromArgb(bytes[i], bytes[i + 1], bytes[i + 2])))
				{
					bytes[i] = bytes[i + 1] = bytes[i + 2] = 0;
				}
			}
			System.Runtime.InteropServices.Marshal.Copy(bytes, 0, bmpData.Scan0, bytesSize);
			image.UnlockBits(bmpData);
		}

		/// <summary>
		/// Excludes given colors from the image.
		/// </summary>
		/// <param name="image"></param>
		/// <param name="colors"></param>
		public static void ExcludeColors(Bitmap image, List<Color> colors)
		{
			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
			int bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			for (int i = 0; i < bytesSize; i += 3)
			{
				var tempColor = Color.FromArgb(bytes[i], bytes[i + 1], bytes[i + 2]);
				if (colors.Contains(tempColor))
				{
					bytes[i] = bytes[i + 1] = bytes[i + 2] = 0;
				}
			}
			System.Runtime.InteropServices.Marshal.Copy(bytes, 0, bmpData.Scan0, bytesSize);
			image.UnlockBits(bmpData);
		}
		#endregion

		#region - Point filter. -
		/// <summary>
		/// Performs image point filtering using prescribed set of pixel coordinates.
		/// </summary>
		public class PointFilter
		{
			public static Bitmap Filter (Bitmap src, List<Point> pts)
			{
				var res = new Bitmap(src);
				var rect = new Rectangle(0, 0, res.Size.Width, res.Size.Height);
				var bmpData = res.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

				var bytesSize = res.Height * bmpData.Stride;
				var bytes = new byte[bytesSize];
				System.Runtime.InteropServices.Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

				var stride = bmpData.Stride;
				var newBytes = new byte[bytesSize];
				foreach (var pt in pts)
				{
					// blacken pixel at point [pt].
					var idx = pt.Y*stride + pt.X*3;
					newBytes[idx] = bytes[idx];
					newBytes[idx + 1] = bytes[idx + 1];
					newBytes[idx + 2] = bytes[idx + 2];
				}

				System.Runtime.InteropServices.Marshal.Copy(newBytes, 0, bmpData.Scan0, bytesSize);
				res.UnlockBits(bmpData);

				return res;
			}
		}
		#endregion
	}
}
